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Abstract We present a hydro-geomechanical model for 
subsurface methane hydrate systems. Our model con¬ 
siders kinetic hydrate phase change and non-isothermal, 
multi-phase, multi-component flow in elastically deform¬ 
ing soils. The model accounts for the effects of hy¬ 
drate phase change and pore pressure changes on the 
mechanical properties of the soil. It also accounts for 
the effect of soil deformation on the fluid-solid interac¬ 
tion properties relevant to reaction and transport pro¬ 
cesses (e.g., permeability, capillary pressure, reaction 
surface area). We discuss a ’cause-effect’ based decou¬ 
pling strategy for the model and present our numerical 
discretization and solution scheme. We then proceed 
to identify the important model components and cou¬ 
plings which are most vital for a hydro-geomechanical 
hydrate simulator, namely, 1) dissociation kinetics, 2) 
hydrate phase change coupled with non-isothermal two 
phase two component flow, 3) two phase flow coupled 
with linear elasticity (poroelasticity coupling), and fi- 


We gratefully acknowledge the support for the first author by 
the German Research Foundation (DFG), through project no. 
WO 671/11-1 

S. Gupta 

Chair for Numerical Mathematics, 

Technical University Munich, 

Boltzmannstrafte 3, 85748 Garching bei Miinchen, Germany 
Tel.: +49-89-28918439 
E-mail: gupta@ma.tum.de 

R. Helmig 

Dept, of Hydromechanics and Modelling of Hydrosystems 
University of Stuttgart, 

Pfaffenwaldring 61, 70569 Stuttgart, Germany 
B. Wohlmuth 

Chair for Numerical Mathematics, 

Technical University Munich, 

Boltzmannstrafte 3, 85748 Garching bei Miinchen Germany 


nally 4) hydrate phase change coupled with poroelastic¬ 
ity (kinetics-poroelasticity coupling). To show the ver¬ 
satility of our hydrate model, we numerically simulate 
test problems where, for each problem, we methodically 
isolate one out of the four aforementioned model com¬ 
ponents or couplings. A special emphasis is laid on the 
kinetics-poroelasticity coupling for which we present a 
test problem where an axially loaded hydrate bearing 
sand sample experiences a spontaneous shift in the hy¬ 
drate stability curve causing the hydrate to melt. For 
this problem we present an analytical solution for pore- 
pressure, which we subsequently use to test the accu¬ 
racy of the numerical scheme. Finally, we present a more 
complex ‘^D example where all the major model compo¬ 
nents are put together to give an idea of the model ca¬ 
pabilities. The setting is based on a subsurface hydrate 
reservoir which is destabilized through depressurization 
using a low pressure gas well. In this example, we sim¬ 
ulate the melting of hydrate, methane gas generation, 
and the resulting ground subsidence and stress build-up 
in the vicinity of the well. 

Keywords Methane hydrate reservoir • hydro¬ 
geomechanical model • kinetics-poroelasticity coupling 
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1 Introduction 

Methane hydrates are formed when water molecules 
form a cage-like structure and trap a large number of 
methane molecules within, forming a crystalline solid 
similar to ice [39]. Methane hydrates are thermody¬ 
namically stable under conditions of low temperature 
and high pressure and occur naturally in permafrost 
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regions or below ocean/sea floors m- If warmed or de¬ 
pressurized, methane hydrates destabilize and dissoci¬ 
ate into water and methane gas. Methane hydrates are 
a very dense source of methane gas. One cubic meter 
of methane hydrate stores approximately 164 standard 
cubic meters of methane gas. Also, the energy content 
of methane occurring in hydrate form is immense, pos¬ 
sibly even exceeding the combined energy content of all 
other conventional fossil fuels [29]. 

The research in methane hydrates stems from a ’three¬ 
fold’ motivation driven by concerns of methane gas ex¬ 
traction and production feasibility, global warming po¬ 
tential and climate change concerns, and inherent geo¬ 
hazards of mining/drilling induced destabilization of 
methane hydrates in subsurface reservoirs. Phenomeno¬ 
logical modeling and numerical simulation of these sys¬ 
tems is thus vital 1) for optimizing recovery techniques 
for extracting methane from hydrate bearing sediments, 
2) for conducting studies and making predictions for 
mitigating bore-hole, local and regional slope stability 
hazards, 3) for sequestering carbon-dioxide in gas hy¬ 
drate, 4) for possible application in natural gas storage 
and transport, and 5) for evaluating the role of gas hy¬ 
drate in global carbon cycle. 

Methane hydrate formations are a fairly complex sub¬ 
surface system characterized by a large number of highly 
interdependent physical phenomena. The typical phys¬ 
ical processes occuring in a stimulated hydrate reser¬ 
voir include 1) hydrate phase change, 2) non-isothermal 
multi-phase, multi-component flow, 3) geomechanical 
deformation of the hydrate bearing sediment, and 4) 
change in the hydraulic as well as the mechanical prop¬ 
erties of the hydrate bearing sediment. Thus, any de¬ 
tailed study of these reservoirs and their possible ap¬ 
plications in energy, environment, and quantification 
of geo-hazards requires the development of 1) multi¬ 
physics mathematical model that accounts for the afore¬ 
mentioned processes and captures their interdependen¬ 
cies, and 2) robust, and efficient numerical tools capable 
of handling multi-physics models and performing large 
scale simulations. 

Several mathematical models have been proposed (e.g. 
Tsypkin [47|, Ahmadi et al. [T], Yousif et al. [50], Sun 
and Mohanty [43], Liu and Flemmings EH, Moridis [SH 
|32]) and different numerical simulators have been de¬ 
veloped (e.g. MH21-HYDRES [52], STOMP-HYD [53], 
UMSICHT HyRes [54], TOUGH-HYDRATE [M]) for 
simulating hydrate reservoirs. These models and simu¬ 
lators consider mainly hydrate phase change and fluid 
flow while neglecting the geomechanical effects. 


It has been widely recognized in the hydrate commu¬ 
nity that the destabilization of hydrates can cause sig¬ 
nificant consolidation and ground deformation, and in 
extreme cases it can even trigger landslides [42]. Thus, 
a lot of experimental work has been done in charac¬ 
terizing the mechanical properties and deformation be¬ 
haviour of methane hydrates and hydrate-bearing sedi¬ 
ments (e.g. Hyodo et al. EHUini, Lee et al. [25]). Several 
mathematical models have been proposed to extend the 
above mentioned hydrate-reservoir model concepts to 
include geo-mechanics (e.g. Rutqvist and Moridis [36] . 
Kimoto et al. [H]). Rutqvist and Moridis [SB] and Klar 
et al. [23] have coupled TOUGH-HYDRATE with the 
commercial geomechanical code ELAC^^ [20] to inves¬ 
tigate the hydro-geomechanical behaviour of hydrate 
reservoirs. Kimoto et al. m have developed their own 
chemo-thermo-hydro-mechanical simulator. Their sim¬ 
ulator uses an elasto viscoplastic model to simulate de¬ 
formation. 

In our work, we focus on building a consistent math¬ 
ematical and numerical framework for hydrate systems 
from Aground zero’ up using a multi-physics approach. 
The primary objective is to capture the dynamic cou¬ 
pling between the transport and mechanical processes 
observed at macroscopic scales. In this paper, we first 
present our mathematical model. The application in fo¬ 
cus is limited to gas production from thermally stimu¬ 
lated or de-pressurized reservoirs. We break the model 
down into its functional building blocks and try to iden¬ 
tify the mechanism of the information exchange be¬ 
tween them. This is important to establish a consis¬ 
tent feedback loop between the processes. It also helps 
us to clearly identify the various couplings present in 
the multi-physics model. We then present the numer¬ 
ical solution scheme followed by test problems where 
we methodically isolate each of the couplings identified 
during the model break-down step. Through these test 
problems we 1) verifying each of the model components 
making up the hydrate simulator, and 2) show the ver¬ 
satility of the model in the variety of hydrate reservoir 
related problems it can handle. Einally, we present a 
3D example problem where a typical subsurface hy¬ 
drate reservoir is destabilized by depressurization. 

Nomenclature 

Xa Mole fraction of component k = H 2 O in 
phase a = g^w 

gCH 4 CH 4 generation rate 

gH 20 generation rate 
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_gHyd^_gh Hydrate consumption rate 
Qh Heat of hydrate phase change 

Volumetric injection rate for phase a = g^w 
€ Volumetric or isotropic strain 

g Gravity vector 

Diffusion flux of component n = CH 2 O m. 
phase a = g^w 

u Sediment displacement vector 

Vg Sediment displacement velocity 

Total phase velocity 

vp Phase velocity relative to the sediment 

ga Dynamic viscocity of phase a = g^w 

Ush Poisson ratio for composite solid 

Actual or total porosity 
(/)e// Apparent or ffective porosity 
Density of phase ^ = g^w^h^s 
psh Density of composite solid 

a , a' Isotropic total and effective stresses 
r Tortuosity 

e strain tensor 

a , a' Total and effective stress tensors 

Ars Specific reaction area 

As Specific surface area of hydrate-bearing 

sediment 

Bulk modulus of phases ^ = g^w^h, s 
Bjn Bulk modulus of bulk REV 
Bsh Bulk modulus of composite matrix 

Cpa Specific heat capacity at constant pressure of 
phase a = g^w 

Cv^ Specific heat capacity at constant volume of 
phase j = g^w^h, s 

Binary diffusion constant in phase a = g^w 
Eh , Es Young’s modulus for hydrate and soil 
Esh Young’s modulus for composite solid 

fg Gas phase fugacity 

Gsh , ^sh Lame’s parameters 
ha Flow enthalpy of phase a = g^w 

H Henry’s constant for methane 

K Intrinsic permeability of hydrate bearing 
sediment 

Thermal conductivity of phase ^ = g^w^h^s 
/Cgjj Lumped thermal conductivity of the REV 
kd Hydrate dissociation rate constant 

kra Relative permeability of phase a = g^w 

kreac Rate coustaut for kinetic phase change of 
hydrate 

Molar mass of component k = H 2 O, CH 4 , Hyd 
Nnydi Hydration number 


Pfj^o Saturated water vapor pressure 
Pa Pressure of phase a = g^w 
Pc Capillary pressure 

Pe 5 Peqb Equilibrium pressure for hydrate phase 

Peff Effective fluid pressure 

Ru Universal gas constant 

Sar Effective aqueous phase saturation 

Sp Saturation of phase ^ = g^w^h 

T Temperature 

Uj Internal energy of phase ^ = g^w^h, s 


2 Mathematical model 


We consider three molecular components: CH 4 , H 2 O 
and CH^. {H 20 )^^ (Hydrate), which are present in three 
distinct phases: gaseous, aqueous, and solid. The gaseous 
phase comprises of molecular methane and molecular 
water in vapour form. The aqueous phase comprises of 
molecular water and dissolved molecular methane. The 
solid phase comprises of pure methane hydrate and soil 
grains. The soil grains are assumed to form a material 
continuum which provides the skeletal structure to the 
porous medium. We shall refer to this as solid-matrix. 
The aqueous, gaseous, and hydrate phases exist in the 
void spaces of this solid matrix (See Fig. [^. At this 
stage the adsorption of methane gas on the surfaces of 
the solid matrix and the hydrate is not considered. 


At the pore-scale, we make a distinction between the 
actual pore-space^ which is the void space outside the 
soil-grains, and the apparent pore-space which is that 
part of the void space that is not occupied by hydrate 
and is thus available for flow of water and gas. At the 
REV-scale, this translates to actual or total porosity 

0 = and apparent or effective porosity (peff = 
Vt 

Vp-Vh 


Vt 


(See Fig. 


1). This distinction is important as it 


gives us a conceptual advantage in isolating the effects 
of deformation and hydrate melting on the hydraulic 
properties of the porous medium. (See Section 2.2.4) 


To describe the hydraulic properties in Section |2.2.4| 
and the effect of hydrate melting on these properties, 
we make an additional assumption that at the pore- 
scale the hydrate coats the soil grain perfectly, and the 
water phase forms a him over the hydrate. 


In the subsequent discussion, the phases occupying 
the pore space (gaseous, aqueous, hydrate) will be de¬ 
noted by ’/5’ = g^w^h respectively, the mobile phases 
(gaseous and aqueous) will be denoted by ’a’ = re, 
and the molecular components will be denoted by ’/i:’ 
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Fig. 1 : Pore-scale to REV-scale 


= CH 4 , H 2 O, Hyd. The soil matrix will be designated 
with the subscript ’ 5 ’. The soil+hydrate composite ma¬ 
trix will be designated with the subscript ^sh\ ’ 7 ’ would 
be used to denote all phases, i.e., ^ h, and s. 

Assumptions The definition of the mathematical model 
is based on the following set of assumptions. 

— Gas hydrate reservoir is of SI type, i.e., the gas caged 
in the hydrate is purely methane. 

— Ice formation is neglected. 

— Aqueous phase is free of any salinity. 

— Local thermal equilibrium (LTE) condition prevails 
in a representative control volume (REV) so that 
the solid and the neighboring fluid are at the same 
temperature. 

— Gaseous and aqueous phases are treated as ideal 
mixtures. It is further assumed that the chemical 
species in the fluid phases attain chemical equilib¬ 
rium instantaneously. It is important to emphasize 
that the chemical equilibrium exists only between 
the fluid phases, not between the fluids and hydrate. 
The hydrate phase change is a relatively slower pro¬ 
cess compared with the methane dissolution and wa¬ 
ter evaporation processes. 

— Elow velocities of mobile gaseous and aqueous phases 
belong to Darcy’s Law regime. 

— Hydrate is of grain-coating type and remains per¬ 
fectly adhered to the pore walls in the soil matrix. 

— Soil and hydrate form a composite material. This 
material (composite-solid-matrix) is treated as a con¬ 
tinuum phase, and the stresses are considered to 
be acting on this composite-solid-matrix as a whole 
(and not on the soil matrix alone). The mechanical 


behaviour of this composite material is described 
using a linear elastic stress-strain constitutive law. 
The mechanical properties of this composite mate¬ 
rial are assumed to vary with the composition and 
stress-state of the composite-solid-matrix. 


2.1 Governing equations 

The transport processes characterizing the gas pro¬ 
duction from a typical sub-surface methane hydrate 
reservoir can be described by invoking the conservation 
laws for mass, momentum, and energy described for the 
macroscale properties of the porous medium m- 

Mass conservation The mass conservation for water 
and gas is written componentwise, i.e., for each compo¬ 
nent n = H 2 O, we get 

^ [dt Pa Xa 5a)] + ^ [V • ((). /)« xS ^a 

a a 

= 5][V-(().5„ JS)] + 5'^ + . (1) 

a a 

The mass conservation for the hydrate phase is given 
by 

dt {4> Ph Sh) + V ■ {(j> ph Sh Vh,t) = 9^ ■ (2) 

The mass conservation for the soil phase is given by 
dt [(1 - 4 >) Ps] + V ■ ((1 - 0) ps Vs) = 0 . (3) 
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Momentum conservation For mobile phases in porous 
medium, under certain simplifying assumptions, the mo¬ 
mentum conservation can be reduced to Darcy’s Law 
|15| . which is stated as, 

V„ = — {VPa -Pag) ■ (4) 

Men 


Here, K is the intrinsic permeability of the composite 
matrix, kr^a and jj^a are the relative permeability and 
the dynamic viscosity of the phase a respectively. 

Vq, is the velocity of the mobile phase relative to the 
soil matrix. The hydrate phase is immobile relative to 
the soil-matrix, i.e., = 0. The total velocity of any 

phase (3 occupying the pores is given by + 

(pSpVg. The soil phase velocity is the rate of deformation 
of the soil matrix, and is given by = dfU. 


Momentum conservation for the composite solid ma¬ 
trix is given by 

V ■ a + psh g = 0 ■ (5) 


Energy conservation For describing the energy con¬ 
servation in the porous medium, one energy balance 
equation is sufficient since local thermal equilibrium has 
been assumed m- The energy balance equation is thus 
given by 


dt (1 - (j)) psUs + ^ ((/) S'/? Pp Up) 

P 

+ Pa Xa Sa ha)] 

a 

= V-KffVT + Qh + ha) 


(6) 


where. 


kiff = {i-4>)kl + Y.Y.^^ 5a K) + 4>Sh K 

a K 




Uj = 



Cpa 


Cvj 


dT 

dT . 


Furthermore, the pressures of the fluid phases are re¬ 
lated through a capillary pressure Pc as Pg — = Pc . 

This pressure difference occurs across the gaseous and 
aqueous phase interface due to balancing of cohesive 
forces within the liquid and the adhesive forces between 
the liquid and soil-matrix. The parametrization used for 
approximating Pc is further elaborated upon in Section 

[2231 

2.2 Constitutive relationships 

In the mathematical model described in Section 12.1! 
following variables can be identified, 

Sp , Xa ■■ Pa , Pc , T , 4> , a , \1 , 3'^ , , Qh 

(7) 

i.e., the total number of variables is 24. (The vectors and 
tensors are considered as single variables and n denotes 
the mobile components, i.e. n = CH 4 ^^ H 2 O.) However, 
the number of governing equations add up to only 12 . 
To close this system 12 additional constitutive relation¬ 
ships are defined for , Pc , d , , and Qh in 

this section. Some other properties which are important 
for modelling hydrate reservoirs are also discussed. 

2.2.1 Vapor-liquid equilibrium 

The two-phase CH^ — H 2 O fluid system is assumed to 
be in a state of vapor-liquid-equilibrium. To calculate 
the concentrations of CH 4 and H 2 O in both gaseous 
and aqueous phases Henry’s Law [3] and Raoult’s Law 
[3] for ideal gas-liquid solutions are invoked: 

For dissolved methane, = H{T) Pg ( 8 ) 

psat fp\ 

For water vapor, = Xw^^ - (9) 

Pg 

Using Eqn. (|^, Eqn. ([^, and the summation condi¬ 
tions for each a, xS = 1 : mole fractions can 

be calculated explicitly. 

In Eqn. (§, H is Henry’s constant for methane gas 
dissolved in water and is calculted in our model us¬ 
ing the empirical relation from NIST standard refer¬ 
ence database m In Eqn. Pfpo saturated 

vapour pressure for water in contact with methane gas. 
We use Antoine’s equation [3] to calculate Pfpo' 


Closure relations The saturations of the phases oc¬ 
cupying the pores satisfy the summation condition 
T,pSp = i ■ 

Additionaly, for each mobile phase a, the constituting 
component mole-fractions also satisfy the summation 
condition XS = 1 • 


2.2.2 Pick’s law for diffusive mass-transfer flux 

The diffusive solute flux through sediment is cal¬ 
culated using Tick’s Law m, as stated below 

JS = -tD^ iPa^Xa) ■ (10) 
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Fig. 2: Methane hydrate P-T phase curve 


The gas-phase binary diffusion coefficient for low den¬ 
sity binary Ci ^4 — H 2 O system is estimated using the 
empirical relationship proposed by Stattery and Bird 
El- For the aqueous phase binary diffusion coefficient, 
Wilke-Chang correlation m for dilute associated liquid 
mixtures is used. 


Ars is the specific reaction area available for the kinetic 
reaction to occur, given by 

= FrAs (14) 


where, Fr is the fraction of the pore surface area that 
is active in hydrate kinetics [44]. The pore surface area 
As is a hydraulic property of the solid matrix (Refer 
Sec. 2.2.4, Eqn. 21). 


The methane gas fugacity fg in the kinetic model is 
computed based on the Peng-Robinson’s thermodynamic 
equation of state for methane [34|. The equilibrium 
pressure for the methane hydrate Pe is determine using 
the Kamath and Holder correlation m as given below, 

Pe = Al exp ^^2 - • (15) 


Hydrate phase change kinetics 

Methane hydrates, upon heating or depressurization, 
decompose to produce methane gas and water as shown 
in Fig. Chemically, this phase change process can be 
expressed as ^ CH 4 FNnyd’ H 20 ^ 

where, N^yd is the hydration number. 


Methane hydrate dissociation reaction is an endother¬ 
mic process. Conversely, methane hydrate re-formation 
is an exothermic reaction. The heat of reaction for hy¬ 
drate phase change is modelled by 

«-£(•-ill) 


This non-equilibrium phase change of methane hy¬ 
drate is modeled by the Kim-Bishnoi kinetic model (8). 
The rate of gas generated or consumed on hydrate de¬ 
composition or reformation is given by 

gCH^ = kreac McHA {Pe “ fg) • (H) 

Correspondingly, the rates of water and methane hy¬ 
drate generated/consumed are given by 

y y McHi 

_gHyd ^ -CHi Mnyd ^ .^ 2 ) 

Mch, ^ ^ 

In Eqn. kreac IS the kiuctic rate constant given 

by 


2 . 2.4 Properties of the fluid-matrix interaction 

Capillary pressure On the macro scale, the capil¬ 
lary pressure in a porous medium is an average pressure 
depending on the pore-size distribution and the aque¬ 
ous phase saturation. Several parameterizations exist 
which relate the capillary pressure and effective aque¬ 
ous phase saturation using soil specihc parameters. Our 
model uses one of the most common parameterizations 
proposed by Brooks and Corey [9]. 

Eor an un-deformed, un-hydrated soil matrix the cap¬ 
illary pressure is expressed as a function of effective (or 
normalized) aqueous phase saturation, as given below 

PcO = Pentry (17) 


kreac = ^reac (13) 

where. Fact is the activation energy, and kreac is the 
intricsic rate constant for the kinetic phase change. 


where, Pentry is the gas entry pressure, \bc is the soil 
specihc parameter depending on the pore-size distribu¬ 
tion, and Swe is the normalized aqueous phase satura¬ 


tion given by S^e = 


{Srur T Pgr) 
Ph {Pwr T Pgr) 
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The effect of presence of hydrate in the soil matrix 
and the changing hydrate saturation on the capillary 
pressure Pc is modelled by scaling ■PcO with a scaling- 
factor which is a function of Sh p^l35] . Also, the 
effect of changing porosity due to deformation of the 
porous matrix is accounted for by scaling Pco with a 
scaling-factor which is a function of (j) using Civan’s 
power-law correlation m- Thus, the capillary pressure 
is given by 


Relative permeabilities The relative permeability 
factors for both mobile phases are evaluated using the 
Brooks-Corey model in conjunction with the Burdine 
theorem m , as. 




2 + 3Xbc 

^BC 


, and 


krg = (1 - U “ 


( 20 ) 




(18) 


where. 


^Pc 

TSh 


(1 - Sh) 


mXBC 


and /|’‘= 


(h { 1 - \ 

(f) yi-cpo) 


Specific surface area An important property of the 
porous matrix is it’s specific surface area, A^, which is 
defined as the ratio of the total internal surface area of 
the pores enclosed within an REV to the total volume 
of the REV. The correlation proposed by Yousif m 
is used for estimating the specific surface area of the 
porous matrix. 


where m and a are model parameters. 


Intrinsic permeability The intrinsic permeability, 
P, is related to the connectivity of the pore spaces and 
the grain size of the soil. It is a property of the soil 
matrix and is independent of the pore-fluids. An esti¬ 
mate of the intrinsic permeability can be made using 
mathematical expressions such as those proposed by 
Bear [6], or Mualem [33]. Usually, however, the intrin¬ 
sic permeability is evaluated experimentally as part of 
the characterization of the soil sample. 


The effect of changing hydrate saturation on the in¬ 
trinsic permeability is modelled by scaling the initial 
or reference intrinsic permeability of the sediment Kq 
with a scaling-factor /|^ which is a function of Sh [121 
l35] . and, the effect of changing porosity due to defor¬ 
mation of the porous matrix is accounted for by scaling 
Kq with a scaling-factor which is a function of (p 
using Civan’s power-law correlation m- Thus, the in¬ 
trinsic permeability for the hydrate sample is modelled 
as. 


K = Ko- {Sh) ■ {4>) 


(19) 


where. 


/£=(l-S.)^ .and = 

These scaling factors (for both. Pc and k) are derived 
based on the assumption that hydrate phase sticks uni¬ 
formly at the pore surface. Eor the ideal case of a spher¬ 
ical pore geometry, m = 3. 


= Y where, ())e// = (/> (1 - S'ft) . (21) 

Hydraulic tortuosity Tortuosity is empirically re¬ 
lated to porosity, as, 

T = (jp where, 1 < n < 3 . (22) 

2.2.5 Poro-elasticity 

Principle of effective stress When a porous fluid- 
filled soil encounters an external load, the stress is partly 
supported by the soil matrix and partly by the pore- 
fluids. The deformation of the porous medium is ef¬ 
fected by only that part of the total stress that is sup¬ 
ported by the soil matrix. This stress, introduced by 
Terzaghi [46], is called the effective stress. Using this 
concept, the total stress a appearing in Eqn. <© can be 
decomposed as. 


^ ^ T ^biotPeffl 


(23) 


where, a is the total stress acting on the bulk porous 
medium, d' is the effective stress acting on the com¬ 
posite skeleton, and Pg// is the effective pore-pressure 
exerted by the mobile phases, given by 


e// 




-Pw + 


:_ p 

Sp 


ai)iot is Biot’s parameter. One of the generally ac¬ 
cepted expressions for aijiot in rock-mechanics applica- 

tions is aijiot = 1 — [Z]. Here, Bsh is the bulk 

Bsh 

modulus of the composite matrix, and Bm is the bulk 
modulus of the porous medium. 
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Stress-strain behaviour The deviatoric stress vs. 
axial strain response obtained from tri-axial tests of 
distributed hydrate-bearing sands reported by Masui et 
al. [28] and Miyazaki et al. m show that the hydrate- 
sand specimens are not elastic bodies, but it is possible 
to consider the stress-strain relationship to be elastic if 
the range of application is sufficiently limited to small- 
strain cases far away from the critical state. 


In our model, we use the linear-elastic constitutive law 
to describe the stress-strain response of the hydrate-soil 
composite matrix, given as. 


Compressibility The grains of the composite mate¬ 
rial are assumed to be incompressible, but the bulk 
material as a whole is compressible. This compressibil¬ 
ity can be attributed to the fact that due to the pore 
pressure variations or isotropic loading the grains re¬ 
distribute, which, on macro scale, manifests as change 
in density of the solid material. This change in density 
can be modelled as. 


d _ Psh 

dt Gsh (1 - 0e//) 





where, a is the isotropic stress. 


(27) 


a' = 2 Gsh e + ^sh{tr e) I (24) 

where, Gsh and Xsh are the Lame’s parameters for the 
elastic composite-matrix and can be defined in terms of 
the apparent elastic mechanical properties (e.g. Young’s 
modulus Eshi and Poisson’s ratio Ugh) as, 

^ _ ^sh A _ ^sh ^sh 

~ 2{l + nsh) ’ + • 

(25) 

e is the linearized strain, given by e = - (Vu + V^u). 

Elastic properties From the tri-axial tests it is observed 
that the presence of methane hydrate, in general, in¬ 
creases stiffness and leads to higher strength. Also, the 
effect of methane hydrate saturation on the Poisson ra¬ 
tio appears to be small. Further details on the general 
trends of mechanical properties of methane hydrates 
can be found in Waite et al. EH, and Soga et al. m- 

To make the model consistent with these observations, 
we assume the Poisson ratio Ugh to be a constant, and 
we define the Young’s modulous using the expression 
proposed by Santamarina and Ruppel (38], given as, 

Esh = Eso(E\ +cEh {Shf (26) 

where, Ego is the isothermal Young’s modulus of hydrate- 
free sand at the reference confining stress of ctco = 1 
kPa, b is the sensitivity of the Young’s modulus of 
hydrate-free sand to confining stress dc, c is the con¬ 
tribution of the isothermal Young’s modulus of hydrate 
for a given pore habit, i.e., pore filling, cemented (grain¬ 
coating) , or patchy, and d is the nonlinear effect of hy¬ 
drate saturation. 


3 Numerical solution 

The mathematical model describing the hydromechan¬ 
ical processes in a hydrate system actually contains 
two sub-classes of models, the flow and transport model 
comprising of the mass, momentum, and energy balance 
equations for the phases occupying the pore spaces in 
the hydrate formation, i.e., equations Q, ([^, Q, and 
([^, and the geomechanical models comprising of mo¬ 
mentum balance equation The soil phase mass bal¬ 
ance, Eqn can be seen as a glue between these two 
sub-models (Fig. [^. 

The interaction between these two models manifests 
physically in the form of, a) changes in the hydraulic 
properties (total porosity and permeability) due to de¬ 
formation of the solid matrix structure, b) shift in the 
seepage velocity of the pore fluids due to the rate of de¬ 
formation of the solid matrix, and, c) changes in the me¬ 
chanical properties of the solid matrix (strength, bulk 
modulus, density, etc.) resulting from the flow dynam¬ 
ics of the pore-fluids and the loss in cementation due 
to melting of the hydrate phase. In other words, each 
model affetcs the other model by altering it’s proper¬ 
ties. Thus, the nature of the coupling between these 
two models is dynamic, but weak. 

We use this observation to our advantage to devise a 
decoupled iterative solution strategy. Broadly speaking, 
we first decouple the flow model and the geomechani¬ 
cal model, solve them separately for a given time-step, 
and then iteratively re-introduce the coupling through 
a block Gauss-Seidel solution scheme. 

We chose the following set of variables as the pri¬ 
mary variables: the gas phase pressure P^, the aque¬ 
ous phase saturation the hydrate phase saturation 
Sh^ the temperature T, the total-porosity 0, and, the 
composite-matrix displacement u. All other variables 
can be derived (explicitly or implicitly) from this set of 
variables. 
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Fig. 3: ’Cause-effect’ based interaction between the flow and geomechanical models. 

The bold arrows represent the direct forward coupling, (e.g. F\ —F\ or Fi — Fh, F 2 —F 3 , F 3 — Fi). The dashed arrows represent 
indirect coupling (e.g. Fi — F 3 ) or backward coupling (e.g. F 2 — Fi, F 3 — F 2 ) 


The system of PDEs comprising the flow-system are 
solved fully implicitly for the variables Pg^ Sh, and 
T. The spatial discretization is done using the fully up- 
winded classical cell centered finite volume method. Or¬ 
thogonal grids aligned with the principal axes are de¬ 
fined and a control-volume formulation with two-point 
flux approximation (TPFA) is used. For time-stepping 
an implicit Euler method is used. 

The geomechanical system is solved for the primary 
variable u. The soil momentum balance equation com¬ 
prising the geomechanical model is spatially discretized 
using the Galerkin finite element (FEM) scheme de¬ 
fined on Ql elements. The FEM formulation described 
by Lewis and Schrefler m is used. 

The soil-phase mass-balance equation is solved sep¬ 
arately for 0. It is spatially discretized using the cell 
centered finite-volume method, and is marched forward 
in time using the implicit Euler method. 

The discretized model, can be represented as a system 
of algebraic equations as 


Fi : Ai (X"+\X")X^+i - Bi(X”+\X") = 0 
which comes from the flow-model. 

Fa ; A 2 (X"+\X")X^+i - B2(X”+\X") = 0 
which comes from the geomechanical-model, and 
F3 ; A3 (X"+\X")X^+i - B3(X”+\X") = 0 
which comes from the total-porosity equation. 

X is the solution vector given as X = [Xi X 2 Xs]^ 
where, Xi = [Pg S^, Sh Tf, Xa = [u]^, and X 3 = 
[0]^. The indices n and n+1 denote the solution at time 
t'^ and respectively. The strong non-linearities in 
each of the sub-systems Fi, F 2 , and F 3 are dealt with 
using a damped Newton-Raphson linearization. Each 
of the resulting linear sub-systems are solved using the 
SuperLU linear solver El This forms the inner iter¬ 
ative loop which takes care of the decoupled solution. 
The coupling between Fi, F 2 , and F 3 is re-introduced 
through an outer iterative loop using a Gauss-Seidel 
scheme, as shown in Fig. 
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initialize iteration counter : k = 0 

_t_ 





r 

solve Fi ^ [JCr'’''""' 



r 

solve Fj 

+ 1 , k+1 ^ 5 ^ 


r 

solve F 3 

j^n + l,k+l j^n+l,k + ljT 


n+ 1, k +1 +1 > k 


set: 

■^n + l , h + l , k + 1 

► A = A 

k^k+1 


r^n vn + 1 


set: X"= X" 
update time : = 


Fig. 4; Block Gauss-Seidel outer iterative solution loop 


problem focuses on the coupling between hydrate phase 
change and fluid flow in the three-phase hydrate model. 
In Section [43l we simulate the classical ID consolida¬ 
tion problem to ensure a correct implementation of the 
poroelastic coupling in our numerical scheme. Following 
this, as an extension to the ID consolidation problem, 
in Section we present a ID test setting where a hy¬ 
drate sample is depressurized while being subjected to 
an external vertical loading. This problem focuses on 
the kinetics-poroelastic coupling. Under simplifying as¬ 
sumptions, an analytical solution for the phase-pressure 
evolution is derived for this setting, which is then used 
to verify the numerical scheme. 

In the final example in Section |4.5| we simulate a 3D 
hydrate reservoir which is destabilized through depres¬ 
surization using a low pressure gas well. In this example 
we put together all the important model components 
including hydrate phase change, non-isothermal effects, 
multi-phase multi-component flow, and poroelastic soil 
deformation. 


4.1 Test 1: Dissociation kinetics model 


The numerical scheme is implemented in the C++ 
based DUNE-PDELab framework [5]. The numerical 
code is flexible in it’s dimensionality and is capable 
of solving problems in ID, 2D and 3D. Eurthermore, 
the block structure of the decoupling scheme makes the 
code modular, and gives the flexibility to solve only 
the flow and transport model or the full geo-mechanical 
model, depending on the problem at hand. 

4 Numerical examples 

The mathematical model described in Section [2] con¬ 
sists of three important parts that are specific to the 
gas-production application. These are the methane hy¬ 
drate dissociation kinetics^ the poroelastic coupling^ and 
the poroelastic-kinetics coupling. In this section, test 
problems which focus on each of these parts separately 
are presented. 

In Section |4.1[ we simulate hydrate dissociation in a 
depressurized lab-scale hydrate sample. In this problem 
the geo-mechanical effects are negligible and reaction 
kinetics dominates over fluid flow. Thus this problem 
effectively isolates dissociation kinetics from the other 
processes. Next, in Section |4.2[ we consider a prob¬ 
lem similar to the five-spot problem in a diagonal flow 
setting, containing a melting hydrate block in the do¬ 
main. The geo-mechanical effetcs are negligible. This 


We consider ID and 2D experiments on hydrate dis¬ 
sociation by depressurization by Tang et al. (2007) [45] 
and Yuhu et al. (2009) [I], respectively. 

4.1.1 ID Case 

Experimental set-up A cylindrical, stainless steel 
cell with internal diameter of 38 mm and length of 500 
mm was used as the main pressure vessel. The cell was 
jacketed with an insulating, impermeable layer and was 
immersed in an air-bath. During each experimental run, 
the dry sands were sieved into size range of 300 — 450 
jam and were pushed tightly into the vessel, resulting 
in a sediment with porosity of 33% and a permeability 
of 300 mdarcy. 

The sediment was saturated with distilled water, and 
the methane hydrate was formed in-situ by slowly in¬ 
jecting methane gas at a pressure higher than the equi¬ 
librium pressure. The hydrate was formed in two stages 
to obtain a homogeneous distribution. 

To perform the dissociation experiment, the back pres¬ 
sure regulator was set to a pressure lower than the 
hydrate equilibrium pressure at the working temper¬ 
ature, and the outlet valve was opened quickly. The 
gas released through the outlet valve was continously 
recorded and a eumulative gas-production curve was 
plotted. 
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Table 1; Initial conditions (Tang et al. (2007)) 


ICs 

Run 2 

Run 3 


[ MPa] 

3.535 

3.584 

Ti 

[°C] 

1.54 

2.08 

Pw,i 

[%] 

29.61 

19.25 

Sh,i 

[%1 

21.83 

25.44 

Ki 

mD ] 

300 

300 


[%1 

30.8 

30.8 


Numerical simulation A schematic of the test do¬ 
main for this experiment is shown in Fig.|^ The domain 
is discretized into 100 cells along the X-axis and the sim¬ 
ulation is performed in ID. The geo-mechanical block is 
switched off and only the flow-transport block is solved. 
Gravity is neglected. The depressurization (i.e. back 
pressure regulation) is considered at the left boundary. 
The left boundary also serves as the gas outlet. Two de¬ 
pressurization modes have been considered for this test. 
In the first {Test-ID:Run2)^ the pressure is decreased 
from 3.535 MPa to 0.93 MPa at Ti,ath = 1.54^0. In the 
second {Test-ID:Run3)^ the pressure is decreased from 
3.584 MPa to 1.94 MPa at Ti^ath = 2.08^C. The total 
dissociation process for Run2 and RunS was reported 
by Tang et al. to last 40 and 110 minutes respectively. 
So, the tend for the numerical simulations was chosen 
accordingly. The initial and boundary conditions are 
listed in Table and Table respectively. 

Tang et al. (2007) used TOUGH-Fx/Hydrate to sim¬ 
ulate the experimental data. The value of the intrinsic 
rate constant as reported by Tang et al. was back 
calculated to fit the experimentai data. We have used 
the same method to caiibrate our modei and obtain the 
best vaiue of for each depressurization mode {Run2 
and RunS). The resuiting vaiues for kinetic-parameters 
AEajR and are listed and compared in Table js} 
^Weported^^ refers to the vaiues reported by Tang, and 
''fitted'^ refers to the vaiues obtained from our caicu- 
iations. Aiso, since the parameterization for Peqb was 
not reported by Tang et ai., we have used the standard 
reiationship proposed by Kamath and Hoider [21] for 
pure methane dissoived in distiiied water. 

Results Fig. and Fig. |6b|show the comparison be¬ 
tween cumuiative gas voiume curves obtained experi- 
mentaiiy and numericaiiy for Run2 and RunS, respec- 
tiveiy. Our numericai resuits show a very ciose overaii 
match with the experimentai resuits. 

4.1.2 2D Case 

The experimentai set-up for hydrate formation and 
dissociation processes used in this experiment are very 



Fig. 5: Test setting for ID experiment (Tang et al. (2007)) 


Table 2; Boundary conditions (Tang et al. (2007)) 


at X = 0, t > 0 

Run 2 

Run 3 

Poutlet [MPa] 

0.93 

1.94 

Tbath [°C] 

1.54 

2.08 


ai X = L, t > 0 


Trig = 0 
TTlw = 0 

VT = 0 


Table 3: Kinetic parameters (Tang et al. (2007)) 



AEajR 



m 

r mol 1 
Lm2.Pa-s J 

Run 2 (reported) 

9400 

1.7 X 10^ 

Run 2 (fitted) 

9400 

1.7 X 10^ 

Run 3 (reported) 

9400 

1.4 X 10^ 

Run 3 (fitted) 

9400 

0.8 X 10^ 


similar to the ID experiment by Tang et al. described 
above. The main difference is the sample geometry, which 
is cylindrical in the ID case and square wafer-like in this 
case. Reaction kinetics is essentially only a time depen¬ 
dent process, and the number of spatial dimensions do 
not directly affect the kinetics. However, testing the ki¬ 
netics model in both ID and 2D geometries ensures that 
the spatial coupling between the different model com¬ 
ponents are correctly resolved, and that no spurious 
spatial effects manifest in the simulation of dissociation 
process. 

Experimental set-up The hydrate formation and 
dissociation unit was a stainless steel vessel with length, 
width, and thickness of 380 mm, 380 mm, and 18 mm 
respectively, and was immersed in an air-bath. The pro¬ 
cedure for sand sample preparation and in-situ hydrate 
formation were similar to that described in Sec. 14.1.11 
The resulting sediment porosity and permeability were 
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Table 4; Kinetic parameters (Yuhu et al. (2009)) 


AEajR 

m 

9752.73 


r mol 1 
\.m^-Pa-s\ 

3.6 X 10"^ 

Qh 

[A] 

-1050 T +3527000 


fRUSSURIZATION 
BOUNDARY LW 


- 38 ci? > 


COMPUTATIONAL DOMAIN 


I Hydrate Sample 

"'Adlabati'cTi'mpermeabr 

NO FLOW BOUNDARY 


Height = 1.8 cm 


P gas=Back 

Pressure 
= 2.25 MPa 
P ~ Pbath 
= 1.7°C 


adiabatic, impermeable boundary 


Sw, =0.389 

Sh, =0.176 K,= 1.97Darcy 

Pg,=3.24MPa 4>i=0.4 

r, =i.7°c 


irig=ni^=0 

P ~Pbath 

= 1.7°C 


adiabatic, impermeable boundary 
Fig. 7: Test setting for 2D experiment (Yuhu et al.(2009)) 


are prescribed at the left boundary. The reaction-kinetics 
parameters are listed in Table 


Results Fig. and Fig. show the experimental 
and numerical comparisons of the gas generation rate 
and the gas pressure for this setting. Our numerical re¬ 
sults show good agreement with the experimental val¬ 
ues, especially towards steady state. 


4.2 Test 2: Three-phase hydrate model 

In this section, we test the coupling between the ki¬ 
netics model and the two phase flow model in our nu¬ 
merical scheme. 


40% and 1.97 Darcy, respectively, and the hydrate sat¬ 
uration was 17.6%. 

To perform the dissociation experiment, the back pres¬ 
sure was reduced from an initial pressure of 3.24 MPa 
to 2.25 MPa, and the outlet valve was opened quickly. 
The bath temperature was maintained at 1.7^C. The 
gas released through the outlet valve was continously 
recorded and a gas production rate was plotted. 

Numerical simulation A schematic of the test do¬ 
main for this experimental set up along with the ini¬ 
tial and boundary conditions are shown in Fig. The 
domain is discretized in 10 x 100 cells, and the simula¬ 
tions are performed in 2D. The geo-mechanical block is 
switched off and only the flow-transport block is solved. 
Gravity is neglected. Depressurization and gas outlet 


For this, we use an artificial setting similar to the five- 
spot problem, with the addition of a gas source in the 
domain in the form of a dissociating block of hydrate. 
This test ensures a correct implementation of the con¬ 
vection, diffusion, and the reaction terms in the 2D 
numerical scheme. It also ensures that the numerical 
scheme does not produce any spurious grid-based ef¬ 
fects. 

Problem set-up A schematic of the test domain is 
shown in Fig. The domain is a unit square with a 
0.3m X 0.3m hydrate block located in the center. The 
domain is initially saturated with water. The hydrate 
saturation in the block is 50%. Point A at (0,0) is the 
gas well. Neumann water-outflux B.C. is prescribed at 
A. The diagonally opposite point B at (1,1) is held 
constant at initial pressure. The rest of the domain 
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(a) Gas pressure profile 


boundaries are closed and adiabatic. The depressur¬ 
ization caused by water outflow at A is expected to 
destabilize the hydrate block causing it to dissociate. 
The released gas must then get drawn towards the low 
pressure in the gas well at point A. 

Numerical simulation For this problem, the ge¬ 
omechanical block is switched off and only the flow- 
transport block is solved. The initial and the boundary 
conditions for the problem are specified in Fig. The 
hydraulic properties, hydrate stability curve parame¬ 
ters, and the dissociation kinetics parameters used in 
the numerical simulation are listed in Table HI The end- 
time for this problem is chosen as tend = 500 minutes. 
For the base test (runO), the domain is discretized uni¬ 
formly into 20 X 20 cells. The time-step size is kept con¬ 
stant at 120 seconds. To check the mesh dependency of 
the numerical scheme, the mesh is successively halved, 

i.e. — 2 i^^)run 0 ’ (^^)run 2 ~ 2 (^^)runl 

, and ^ • The time-step size is 

also successively halved so that AxjAt ratio remains 
constant for each of the test-runs. 


8000 



Time (min) 

(b) Gas generation rate 

Fig. 8: Results (comparison with Yuhu et al. (2009)) 


Results The gas plume takes about 300 minutes to 
reach the gas well at A. Fig. pT| shows the screenshots 
of gas saturation in the domain at 100, 200, and 300 
minutes. Fig. 11 shows the line-plot of Sg at 200 and 300 
minutes along the diagonal aligned with flow direction, 
i.e. line X — Y = 0. 


In Fig.pT a ]the solution shows convergence with mesh- 
rehnement. The flow in the right half of the domain (i.e. 
T + X — 1 > 0) is diffusion dominated, whereas, that in 
left half (i.e. T + X — 1 < 0) is convection dominated 
being strongly influenced by the low pressure in the gas 
well. The gas front is more diffusive on a coarse mesh, 
but gets sharper as rehnement is increased. Fig. |llb| 
shows the saturation of the gas plume that reaches the 
gas well at 300 minutes. 


We would like to point out that in Fig. |lla| and Fig. 
fiibl what appears to be a kink in gas saturation at the 
corner of the hydrate zone is not a numerical artifact. 
This kink is caused because of the following physical ef¬ 
fects: the gas velocity in the hydrate free zone is higher 
than that in the hydrate zone (due to difference of al¬ 
most an order of magnitude in the permeabilities) (See 
Fig. 12). This causes the gas to be sucked out of the 
hydrate zone faster than the time required by gas to 
equilibriate inside the hydrate zone. So, the gas begins 
to deplete along the edges of the hydrate zone. Further, 
the extent of the depletion is higher where is lower 
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0.003574 



0.003574 



(a) Runl 

[Ax = , At = 60s) 


(b) Runl 

[Ax = , At = 60s) 


(c) Runl 

[Ax = , At = 60s) 



0.003574 


(d) Run3 

[Ax = , At = 15s) 

Time = 100 min 



0.003574 


(e) Run3 

[Ax = , At = 15s) 

Time — 200 min 
Fig. 10: Five-spot test - Sg profile 



0.003574 


(f) Run3 

[Ax = Y^™ 5 = 15s) 

Time = 300 min 


P= 2.85 MPa 


o 

o 


a 

c- 

a 


no flow, adiabatic 

-1 Pf = 2.85 MPa Ki = 1X10 
:0 T; = 2.25°C cp = 0.3 

_ ^(0.65,0.65) 

= 0.5 

S, ’,i = 0.5 
P^ = 2.85 MPa 

T. = 2.25°C 
K,. = 2.27X10“^", 

(p = 0.3 


B ( 1 , 1 ) 


a 

o 

o 


A ( 0 , 0 ) 


a (0.35,0.35) 


no flow, adiabatic 


m 


Table 5: Five-spot test for hydrate model 

Model parameters 


Brooks-Corey parameters 

Pentry 

[Pa] 

5000 


^BC 

[°C] 

1.5 


Hydrate stability curve 

Peqb 

[kPa 

1 exp ^38.980 

if T > 273.15 



exp (14.717 

if T < 273.15 

Dissociation kinetics parameters 

AEalR 


9400 



r mol 1 Q n N/ 1 n4 


L • 

Pa-si ^ 



Fig. 9: Test setting for 5-Spot test for hydrate model 


(i.e., where Pg is higher). This effect is more clearly vis¬ 
ible in Fig.|13a|and Fig. |13b| which show the Sg profiles 


along X-axis (at F = 0.5 m) at times t = 300 min and 
t = 500 min respectively. The pressure in the right half 
along the X-axis is higher than that in the left half, 
causing to be lower in the right half. Therefore, 

the extent of gas depletion is higher in the right half of 
the hydrate zone. 
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Legend: 

- ref=l, - ref=2, - ref=4, - ref=8 

Fig. 11: Five-spot test for hydrate model 
Sg profile along diagonal line X = Y 

4.3 Test 3: Poroelastic coupling 

In this example we ignore the methane hydrates in 
the medium, thus reducing the model to a simple two- 
phase hydro-mechanical system. We consider the clas¬ 
sical ID consolidation problem by Terzaghi [46] to test 
the fluid pressure response generated by the mechanical 
compression of the soil. This test was originally formu¬ 
lated by Terzaghi for analyzing the time delay observed 
when compressing clay layers and is now considered as a 
standard benchmark test for the coupling relationships 
between fluid and mechanical systems. 

Problem statement The problem set-up consists of 
a confined soil sample surrounded by a circular ring and 
placed in a container filled with water. The sample is 
loaded by a constant or ramped vertical stress at its up¬ 
per surface, and the deformation is measured. The lower 
boundary is impermeable, and the upper boundary is 
fully drained. This is called a confined compression test 
or an oedometer test. Fig. [T^ shows a schematic for this 



Fig. 12: Five-spot test - profile (at t = 30 minutes) 
Red vectors represent in the hydrate zone 
Black vectors represent Vg in hydrate free zone. 


Table 6: Solid-matrix properties for Terzaghi problem 


Property 

Symbol 

Unit 

Value 

Drained bulk 
modulus 

Bs 

GPa 

8.0 

Poisson ratio 

V 

- 

0.20 

Porosity 


- 

0.19 

Permeability 

K 


1.9 X 10-13 

Biot constant 

a 

- 

0.8 


Table 7: Fluid properties for Terzaghi problem 


Property Symbol Unit Value 


Wetting fluid 


Bulk modulus 

Bw 

GPa 

2.933 

Density 

Pw 

kg 

997.05 

Viscosity 

Pw 

Pa • s 

8.9008 X 10 

Non-wetting fluid 
Bulk modulus 

Bnw 

GPa 

1.187 

Density 

Pnw 

kg 

m3 

997.05 

Viscosity 

P^nw 

Pa • s 

8.9008 X 10 


problem. It is expected that the compression of a soil 
sample will be accompanied by an expulsion of pore flu¬ 
ids from the sample. Also, if the soil permeability is low, 
this may take considerable time. In Terzaghi’s original 
work, the pore fluid and the soil particles were both 
assumed to be incompressible, so that the only mech¬ 
anism of deformation was a rearrangement of the par¬ 
ticles. However, Biot’s more generalized consolidation 
framework, which is also the basis of our poroelasticity 
model, accounts for both fluid and soil compressibili¬ 
ties. So, in the further discussion, the fluid and the soil 
are treated as compressible mediums. 
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(a) Time = 300 minutes 



(b) Time = 500 minutes 





Fig. 14: Schemmatic of Terzaghi problem 


Fig. 13: Five-spot test 
Sg profile along X-axis (at Y = 0.5 m) 


The mathematical description of such a problem in ID 
reduces to a fluid diffusion equation of hydrogeology, 


d,P-c 


dz‘^ 


(28) 


where, P is the mean fluid pressure given by P = 
SwPw + SnwPnwi ^ud c is 1-D fluid diffusivity. 


^ = (2m + l)7r/ (21/), L is the total column length, 
and 2 ) is the location in the column downward from the 
applied stress. Hy is the 1-D Skempton coefficient given 
by, 


Hy 


6P 


S(T^, 


a 


^xx — ^yy —C — 0 


H sv Py 


(31) 


where, Bsv is the uniaxial drained bulk modulus, and 
Sy is the 1-D specific storage given by Sy = KK/j^c). 

..... 1 .... . 1 1 A 1 1 ^ 

/i IS the fluid mobility given as — = - I-^-I 

fJj Z \IZyj l^nw / 


For a vertical load ramped linearly at the top 
boundary at a rate dazzidt = &z^ the analytical so¬ 
lution for the pore pressure response is given as 


The complete derivation of the analytical solution (Eqn. 


29) can be found in many of the textbooks on soil me¬ 


chanics, for example, Verrujit (2013) 




32 


L-z 

L 

_ 




cos [//> (L - z)] I (29) 


where, Pq is the total pressure generation given as, 

Po = ^ (P„a,) . (30) 

zc 


Numerical simulation For benchmarking, we use 
the test setting described by Kolditz et al. [24]. A soil 
column of 50 m is chosen. The properties of the rock 
material are listed in Table and that of the two fluid 
phases are given in Table The column is discretized 
uniformly into 200 grid cells. The initial fluid pressure 
in the column is null, and the initial fluid saturations 
are Syj = 0.8 and Snw = 0.2. The hydraulic proper¬ 
ties are chosen as Pc = 0 and ky^yj = ky^nw = 0.5. As 
shown in Fig.[^ a load cFzz = 10 MPa is applied at the 
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top boundary at a loading rate of = 0.0land 0.001 
MPa/s. The bottom boundary of the column is sub¬ 
jected to roller displacement BC and the top boundary 
is allowed to compress freely under the applied load 
The top boundary is a free-drainage boundary. All 
other boundaries are no-flow. 

Results The results of the numerical simulation for 
the two loading rates (d^ = 0.01, 0.001 MPa/s) are 
presented in Fig. |15a| and Fig. |15b[ Compression of the 
column leads to a rapid pressure increase followed by 
dissipiation of pressure over time from top of the col¬ 
umn. It can also be observed that for a lower loading 
rate, the pore-pressure equilibrates faster, whereas, for 
a higher loading rate the pore-pressure takes longer to 
dissipate. The numerical results show a very good agree¬ 
ment with the analytical solution. 


4.4 Test 4: Kinetics-poroelastic coupling (KPE test) 


We now extend Terzaghi’s ID consolidation problem 
to include hydrate kinetics in the poroelastic coupling. 
We consider a confined soil sample which is uniformly 
hydrated and fully saturated with water. A constant 
vertical stress is applied at the top boundary while the 
lower boundary is held fixed. The upper boundary is 
fully drained, while at the lower boundary the initial 
pressure is maintained at all times. For time t < 0“, the 
thermodynamic state of the sample lies on the hydrate 
stability curve, so that P^~ = P ^~, and the hydrate in 


the sample is stable (see Fig. 16). Here, P indicates the 
phase pressure, and Pe indicates the equilibrium pres¬ 
sure for hydrate stability. At time t = 0, the hydrate 
stability curve experiences an instantaneous shift such 
that PO > PO- (while PO = po-). The hydrate be¬ 
comes unstable and begins to dissociate. This generates 
excess pore-pressure which prevents the full consolida¬ 
tion of the sample. 


The schematic for this problem is shown in Fig. pT| 
Although highly simplified, this problem helps us to 
isolate the poroelastic-kinetic coupling^ thus providing a 
framework for validating the numerical implementation 
of our hydro-mechanical code for the hydrate reservoir 
model. 


4-4'^ Prohlem statement 

For this problem we make the following additional as¬ 
sumptions: 

— Gas does not dissolve in water phase, and water 
vapor is not formed, i.e., — 




numerical (^ 
numerical (|^ 
numerical (|^ 
numerical (|^ 


Legend: 

= 0.60) analytical (|; = 0.60) 

= 0.40) analytical 0.40) 

= 0.20) analytical 0.20) 

= 0.04) analytical 0.04) 


Fig. 15: Terzaghi benchmark test results 
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Fig. 16: KPE test - Hydrate stability curve shift at t = 0 



Fig. 17: KPE test - Problem schematic 


1. Based on this assumption, we rewrite the mass 
balance equations for water and methane in Section 
phase-wise instead of component-wise. 

— All the processes, including hydrate phase change, 
are isothermal. 

— There is no suction pressure between the two mo¬ 
bile phases, so that Pg — = 0. Since the phase 

pressures are now equal, we drop the subscript and 
assign the symbol P to the phase pressures through¬ 
out this section. 

— Relative permeabilities are kr^g = kr^w =0.5. 

— Effect of gravity is neglected. 


— Effect of porosity and hydrate saturation on intrin¬ 
sic permeability K is neglected, i.e. K is constant. 


Eurther, we simplify the hydrate reaction kinetics model 
as 

Qg = ko As Mg {Pe - P) 

A, = ^,,0 Sh (32) 

where ko is the rate of hydrate dissociation, and 
is the specific surface area of the hydrate free sample. 
Both ko and Ag^o are assumed to be constant. 


Using these assumptions, we can reduce the mathe¬ 
matical model described in Section [2] to an ODE for 
the pressure P, given by 

+ ^ {Pe - P) ■ (33) 

at at fjjf 

A detailed derivation is given in [Appendix A 


The term S, called the Storativity, is given as 


5 = 




+ 


Ba 


+ 


( a - (pe ) 
Bsh 


e n I AT ^^^^9 . 

term C = Nh -1- - - ko Agfi , 

V Pw Pg Psh J 

Cy is the volumetric strain given as = V • u , and 

/if is the fluid mobility given as — = x (-^- 

Pf ^ xPg Pw 

a is Biot’s constant. 


Eqn. (33) is the storage equation. In this form it can 
be interpreted as: on the REV scale, the compression of 
the soil consists of compression of the pore-fluids and 
the compression of the solid particles, plus the total 
volume of fluid expelled from the REV and the fluid 
generated in the REV. 


Eurther, in Eqn. (33), the term 


ehanical part, and the term 


part. The terms 


dt 


d 

V - — VP 
k-f 


is the me- 
is the flow 


and [C Sh (Pe — P)] are the 


eoupling terms, the former for the coupling between 
the flow and the mechanical models, and the latter for 
the coupling between the flow and the reaction kinetics 
models. 


Eor the ID problem under consideration, we can rewrite 
Eqn. (33) as 


aP, + S^P=— Ap + c Sh{Pe-P) 

dt dt jif dz^ 


( 34 ) 
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In case of ID deformation, the volumetric strain equals 
the vertical strain and is induced by the vertical stress 


d _ d , _ ( d d 

dt"^~ “dt 


(35) 


where, Cm is the compressibility of bulk porous mate¬ 
rial, such that Cm = • 

Bm, 


d 


Thus, we eliminate — in Eqn. (34) using Eqn. (35), 

dt — — 

which gives 
d -p OiCm d 

Jt a^Cm + Sdi'^" 


(f 


+ 


c Sh 


a'^Cra + 5 


+ dz^ 

{Pe - P) . (36) 


At time t = 0 , an external load q is instantaneously ap¬ 
plied, and the equlibrium pressure of hydrates is instan¬ 
taneously changed from P^~ to P^. Since both these 
processes are instantaneous, no fluid is mobilized at 
_ (f 

t = 0, i.e., in Eqn. (36), -r^P = 0. So, from Eqn. 
— dz^ 

(36), we get the initial condition of the sample as 

(^Cm 


t = 0 : p = p^ = 


a^CmPSpCSh 
CSh po 
a-^Cm + S + CSh 


q 

pU 


(37) 


Eor t > 0, the external load remains constant, so 

— azz = 0. The equilibrium pressure also remains con- 
dt 

stant, i.e., Pg+ = P^ = Pq. Thus, from Eqn. (36), 

d ^ K d? ^ 

f z> 0 ' —P =- - P 

dt jjLf {a^CmS) dz^ 


a 

c Sh 

a^Cm + S 


(Pe - P) 


t > 0 


Cr 

|p = C„^F + ft(F.-P) 


(38) 


Cy is the consolidation parameter which comes from 
the Terzaghi’s classical theory of consolidation. Cy is 
the reaetion parameter. It is indicative of the damping 
of the normal consolidation due to dissociation kinetics. 

The boundary conditions at the top and bottom of the 
sample are 


t > 0, z = L : -—P = 0 
dz 

t>0, z = 0 : P = P° . 


(39) 


The ODE in Eqn. (38) is a non-homogeneous ODE. We 
can homogenize it by choosing a new primary variable 
P such that P = Pe — P. Then the initial-boundary- 
value problem (IBVP) can be formally summarized as 


0 < 2 ;<I/, t >0 


— P — C ^ P — C 

dt-^^dz^^ 


z = 0, t > 0 
z = t > 0 
0<z<L^t = 0 


P = Pe- P^ 
0 

P = P^- P^ 


tP 


where, 

0 aCg, q + Pe 

+ S0 + CS° 


(40) 


which is a homogeneous ODE with non-homogeneous 
boundary conditions. An analytical solution can be ob¬ 
tained for this problem using any of the standard tech¬ 
niques for solving ODEs. The final solution for P can 
be written as 


Pe—P{z,t) cosh {0 {L — z)) 
Pe-P^ ^ cosh (OL) 






exp [-Cy {xl + d'^)t]^ 


sin (Xnz) 


(41) 


4.4-2 Numerieal simulation 


Test setting We chose a sample of length L = 1 m 
containing 30% hydrate by volume. The sample is ini¬ 
tially fully water saturated and is contained in a pres¬ 
sure vessel at = 6 MPa. A constant external load 
g = 10 MPa is applied at the top boundary, i.e., at 
2 ; = P = 1 . At the bottom boundary, i.e., 2 ; = 0 , the 
pressure is held constant at the initial value. 

Eig. shows the domain specifications, the initiai 
conditions, and the boundary conditions. The materiai 
properties are iisted in Tabiej^ 


The storage equation governing this probiem, Eqn. 


(40), contains two parameters: Cy and Cy. To test the 
numericai impiementation, we chose three different vai- 
ues of Cy and each with a different order of magni¬ 
tude. Therefore, we run nine tests with aii combinations 
of the chosen Cy and Cy. We controi the parameter Cy 
by varying the sampie permeabiiity n and the param¬ 
eter Cy by varying the dissociation rate constant k^. 
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0 ^ 2 ( 1 , t)=q = 10 MPa 


z 


=L=lm 

i 


z=0 


LC.s 

^h,0 ~ 

^w,0 ~ 

(|) = 0.3 


0000 

u,(0,t)=0 


dz 


P{L,t)=0 


P{0,t)=P, 


Table 9: KPE test - Control parameters 


test 

ID 

K, 

[mD] 

r mol 1 

a 

Cr 

Pe 

[MPa] 

[ ^2 g J 

1 

0.1 

360 

1.53755 

0.289504 

19.151 

2 

0.01 

360 

0.153755 

0.289504 

19.151 

3 

0.001 

360 

0.0153755 

0.289504 

19.151 

4 

0.1 

3600 

1.53755 

2.89504 

7.315 

5 

0.01 

3600 

0.153755 

2.89504 

7.315 

6 

0.001 

3600 

0.0153755 

2.89504 

7.315 

7 

0.1 

36000 

1.53755 

28.9504 

6.132 

8 

0.01 

36000 

0.153755 

28.9504 

6.132 

9 

0.001 

36000 

0.0153755 

28.9504 

6.132 


*a = exp {-AEa/{RT)), where, AEalR = 9400 K and T = 283.15 K 


erated fluids will mobilize. Conversely, the slower the 
dissociation reaction, the higher margin we get for rais¬ 
ing Pe without instantaneously mobilizing the fluids. 


Fig. 18: KPE test - Computational domain settings 
Table 8: KPE test - Material properties 


Simulation and results The domain is discretized 
into 400 cells in z-direction, and the problem is solved 
in IP. The time-step is kept constant at t = 0.1 s and 
the simulation is run until tend = 60 s. 


Property 

Symbol Unit 

Value 

Water phase 

Density 

Pw 

kg • m~'^ 

997.05 

Molar mass 

Mw 

kg ■ mol~^ 

0.018 

Dynamic 

P'W 

Pa • s 

8.9008 X 10-3 

viscosity 




Bulk modulus 

Bn, 

GPa 

2.933 

Gas phase 

Density 

Pg 

kg • m~’^ 

0.717 

Molar mass 

Mg 

kg • mol~^ 

0.016 

Dynamic 

P^g 

Pa • s 

1.0245 X 10-^ 

viscosity 




Bulk modulus 

Bg 

GPa 

0.1013 

Hydrate phase 

Density 

Ph 

kg • m~'^ 

900 

Molar mass 

Mh 

kg • mol~^ 

0.119 

Hydration 

Nh 

- 

5.75 

number 




Young’s 

Eh 

GPa 

1.35 

modulus 




Soil phase 

Density 

Ps 

kg ■ m ^ 

700 

Surface area 


0 riP 

10^ 

Young’s 

Ee 

GPa 

0.3 

modulus 




Solid composite 

Poisson ratio 

Esh 

- 

0.2 

Biot constant 

a 

- 

0.8 


In Fig. the numerically computed pressure values 
for each test case 1 — 9 are compared with the analytical 
pressure P{z^ t) obtained from Eqn. For each test 
case, the pressure solutions are plotted over time at the 
observation points z = 1 m, 0.8 m, 0.6 m, 0.4 m and 
0.2 m. The plots show a good agreement between the 
numerical and the analytical solutions signifying that 
the poroelastic - reaction kinetics coupling terms are 
correctly handled in the numerical code. 


testID-3 The pressure build-up along the length of 
the sample is plotted for testID 3 in Fig. Also, for 
testID 3 a grid-convergence study is performed. The 
mesh is refined from n = 25 cells up to n = 800 cells, 
and correspondingly, the time-step size is reduced from 

Sz 

St = 2s down to St = 0.0625s such that the ratio 

St 

remains constant. The L2 error is calculated at t = 10s 
for each refinement and is plotted against the number of 
grid-cells n on a log-log graph in Fig. 21 It can be seen 
that the error decays linearly with refinement. This is 
in line with the finite volume discretization technique. 


For each test, the value of Pe is chosen such that the 
initial condition of no-drainage is satisfied. The control 
parameters for each of the nine tests are listed in Table 
in It can be observed in Table [9] that as the reaction rate 
constant ko increases the value of equilibrium pressure 
Pe decreases. This is due to the no-drainage condition 
at t = 0. The faster the dissociation, the more the gen- 


4.5 Test 5: 3D hydrate reservoir problem 

So far we considered examples which focussed on sys¬ 
tematically isolated couplings and model components. 
In this section, we present a more complex example 
where we simulate the hydro-geomechanical processes 
in a subsurface hydrate reservoir which is destabilized 
by depressurization using a low pressure gas well. This 
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Fig. 19; KPE test - Comparison of numerical and analytical solutions 
Note: For TestID 5, pressure profiles at 2 ; = 1 m and 0.8 m are equal, for TestID 6 and TestID 8, pressure profile at 2 ; = 1 m, 
0.8 m, 0.6 m, and 0.4 m are equal, and for TestID 9 pressure profile at 2 : = 1 m, 0.8 m, 0.6 m, 0.4 m, and 0.2 m are equal. 


example puts together all the important components of 
our model including dissociation kinets, non-isothermal 
effects, multi-phase multi-component fluid flow, and geo¬ 
mechanics, and qualitatively shows the effects and counter 
effects of various physical processes occuring in the hy¬ 
drate reservoir. The objective of this example is to give 
a first idea about the capabilities of our hydrate reser¬ 
voir model. A detailed quantitative analysis of the prob¬ 


lem and parameter sensitivity study is however beyond 
the scope of this paper. 

Test setting We consider a scaled down 3T) reservoir 
with dimensions 10m x 10m x 5m, as shown in Fig. 
The hydrate is homogeneously distributed in a 4m thick 
layer lying between 0.5m < 2 : < 4.5m, and has a satura¬ 
tion of 40% by volume. The reservoir is fully saturated 
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Pressure (Pa) 

Fig. 20: KPE test - testID 3: Time-evolution of the numerical 
solution of P along the height of the sample 



Fig. 21: KPE test - L2-error vs n 


Constant pressure 
I (10,10,z) 



{0m,0m,0m) 


Fig. 22: Schemmatic of depressurized 3D-reservoir example 


with water and has an initial pressure of 10 MPa. The 
reservoir is depressurized through a low pressure gas 
well located at (0, 0, z). The pressure in the gas well is 
maintained at Pweii = 4 MPa. A constant vertical load 
of 10 MPa is acting on the top boundary of the reservoir 
(i.e. at z = 10 m). The initial and the boundary condi¬ 
tions are listed in Table pT| and Table [m respectively. 


Table 10: Initial conditions for depressurized 3-D-reservoir 
example 


Hydrate layer 

Peff,i = 10 MPa 
Sh,i = 0.4 

at t = 0, for Sg^i =0 

0 < X, ^ < 10 , 0.5 < ^ < 4.5 Ti = 10 OC 

Ki = 0.0198 mD 
= 0.18 


Hydrate-free layers 

at t = 0, for 

0 < a:, ^ < 10 , 2 ; < 0.5 
and 

0 < a:, ^ < 10 , 2 ; > 4.5 

Peff,i = 10 MPa 
Sh^i = 0 

Sg,i =0 

Ti = 10 OG 

Ki = 0.1 mD 

4^eff,i — 0.3 

Table 11: Boundary conditions for depressurized SD- 
reservoir example 

FLOW model 


Gas well at 

x = 0, ^ = 0, 0<^<5 

Pg = A MPa 
= 0 

V ■T = 0 

Pressure constraint at 

X = 10, ^ = 10, 0 < 2 : < 5 

Peff = Peff,i 

Su) — Sw,i 

T = Ti 

No-flow and adiabatic conditions 
on remaining boundaries, i.e.. 

Vg ‘ h = 0 
v-u; • n = 0 

V-T = 0 

GEOMEGHANIGAL model 

Top boundary 

0<x<10,0<^<10,^ = 5 

Bottom boundary 
0<x<10,0<^<10,^ = 0 

Remaining boundaries 

Gzz — 10 MPa , 

^xy — ^yx — 0 

Uz = 0 , 

^ xy — ^yx — 0 

Ux — '^y — 0 , 

(Jzz = 0 


The material properties and other model parameters 
are listed in Table fT^ 


Numerical simulation and results The domain is 
discretized into 30x30x15 cells. Full hydro-geomechanical 
model is solved. To reduce the computational costs the 
decoupling strategy and iterative solution scheme de¬ 
scribed in Section is used. The primary variables be¬ 
ing solved for are: gas phase pressure P^, aqueous phase 
saturation hydrate phase saturation Sh, tempera¬ 
ture T, total porosity 0, and displacements u. Some of 
the other important secondary variables which are cal¬ 
culated as post process include gas saturation Sg, effec¬ 
tive porosity (peff^ intrinsic permeability K, stresses d', 
strains e, etc. The simulation is run until tend = 24 hrs 
with a time step size of dt = 200 s. Selected profiles 
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Table 12; Material properties and model parameters for de¬ 
pressurized 3D-reservoir example 


Thermal conductivities 


kl -0.886 X 10-2 

+ 0.242 X 10-3 T 

i 

1 

1 

-0.699 X 10“® r2 +0.122 x 10“® T® 


0.3834 ln(T) - 1.581 W-m-^-K-^ 

K 

2.1 

ly.m-i -7^-1 

K 

1.9 

W-m-i-K-i 

Specific heat capacities 

Cpw 

4186 

J-fcg-i -K-i 

Cvw 

Cpw + Rh^o j -kg ^ ■ K 1 

Cvh 

2700 

J-feg-i -7^-1 

Cvs 

800 

J-fcg-i -7^-1 

Dynamic viscosities 


273.16+162\ , T \l-t 

’ Pa ■ s 


T + 162 ) V273.15/ 

Hyj 0.001792 exp [ 1.94 4.802^^i® 

Pa • s 


+6.74(27^15)2] 



Densities 


Pg 

Po 

ZRgT 

kg • ra~^ 

Pw vapour: 0.0022^ 

kg • m~^ 

liquid: 1000 

kg • m~^ 

Ph 900 

[i + (“r-U) Vr 

kg ■ m~^ 


sh + (2/3) Gg/i \ V-u 1 


\ 

PsH J l-0e//J 


Ps 2100 

[i+(“r-U) Vr 

kg • ra~^ 


sh + (2/3) Gg/j \ V-u 1 


\ 

PsH J l-0e//J 


Hydraulic properties 


^BC 

1.2 


Pentry 

50 

kPa 

m,a (Eqn. (p^Sf) 

3, 2 


Hydrate kinetics 


uO 

'^reac 

3.6 X 10^ 

mol • m~‘^ • 
Pa~^ • s~^ 

AEajR 

9752.73 

K 


Nnyd 5.75 

Ai, A 2 , As (Eqn. 1000, 

38.98, 

8533.8 

Hi, B 2 (Eqn. @) 56599, 

16.744 

Fr (Eqn. ([^) cj) Sh 

Poroelasticity parameters 


O^biot 

0.6 


^ sh 

0.2 


EsO 

0.3 

GPa 

Eh 

1.35 

GPa 

b, c, d (Eqn. (|26|)) 

0, 1, 1 



showing the state of the reservoir at tend are shown 
in Fig. Fig. |23a| and Fig. |23b| show the melted hy¬ 
drate and the accumulated gas in the vicinity of the 
gas well. Fig. |23c| shows the decrease in temperature 
due to the endothermic nature of hydrate dissociation. 
Fig. |23d| shows the stress built up in the region aroung 
the well where the hydrate is dissociating. Fig. |23e| and 
Fig. |23e| show the change in effective porosity and in¬ 


trinsic permeability as a result of hydrate dissociation 
and soil deformation. The vectors in Fig. show the 
displacements u. The domain is warped with respect 
to displacement to show the ground subsidence around 
the well clearly. The warping of the domain is achieved 
through post-processing using PARAView [2]. 

5 Concluding remarks 

In this article, we have presented a model concept for 
multi-phase, multi-component flow through deformable 
methane hydrate reservoirs. This forms the core of our 
hydrate numerical model and contains only those model 
components which are necessary to simulate the most 
important hydro-geomechanical processes observed in 
a subsurface hydrate reservoir. The structure of this 
model is such that the core can be modularly extended 
to enhance each of the model components to the de¬ 
sired level of complexity depending on the application 
at hand. 

Our focus so far has been to ensure that the dynam¬ 
ics of the hydro-geomechanical interactions are consis¬ 
tently accounted for in our mathematical model. We 
have identified the important physical processes and the 
cause-effect based couplings. Based on this we have pre¬ 
sented our decoupling strategy, where we have discussed 
how we breakdown our complex multiphysics model 
into relatively simpler sub-models. We have also dis¬ 
cussed our solution strategy, which involves first soiving 
the sub-modeis separateiy to obtain a decoupied soiu- 
tion, and then reintroducing the coupiings iterativeiy. 
Through the numericai exampies, each of which isoiates 
an important physicai process in a hydrate reservoir, we 
have shown that our modei is versatiie and is capabie 
of capturing the important coupiings effectiveiy. 
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Fig. 23: Numerical example - Depressurization of a 3D hydrate reservoir 
Selected profiles at time = 2.5 hours 
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Appendix A : Derivation of storage equation 
(Refer Eqn. (33)) 


We re-write the mass conservation equations for gas 
and water phase-wise, i.e., for each fluid-phase a = g,w, 


0 _l_ “1“ ^ oiNcx,t) j.oc 


(Al) 


where, is the volumetric source term for phase a 
given as (xS 5'")- 


Expanding the partial derivatives in Eqn. and 

rearranging gives 


• Vpa j f~Pa^ (05'a) 


P pa^ • = Qa 

4^^a~^Pa Pa~^ {4^^a) Pa^ ' ~ Qa • 


(A2) 


The rate of change of the fluid density is defined as 

—pa = ^^-rPa : where, Ba is the fluid-phase bulk 
at at 


modulus. Using this definition in Eqn. (A2) and divid¬ 
ing by Pa, we get 


^ ^ + V • ^ . (A3) 

Ba at ot pa 


Since we assume Pc = 0, the phase pressures are equal. 
So we drop the subscript and assign the symbol P to 
the phase pressures. 


Next, we sum Eqn. (A3) over a 


g,w, which gives 





+ V ■ + V ■ {(f>SgVgP 


Qw 

Pw 


k 

Pg 


(A4) 


We define the effective fluid-phase saturation Sa,e as 
the volume of fluid phase a = g,w m. the effective pore 
space which is characterized by the effective porosity 


(fe- So, Sa,e — 




, and (fe = 0(1 — Sh) . Also, 


by the summation relationship, Sw P Sg = 1 — Sy , 


or 


1-P. 1-Pi 


= 1 . 
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Substituting these definitions in Eqn. (A4), we get Eqn. (A9) expresses the total mass balance for the 


Eqn. (A5) which we will call the fluid mass balance 
equation. 


whole porous medium consisting of phases j s. 

We now substitute the simplified constitutive relation- 




ships listed in Table [T^ in Eqn. (|A9|), 


+ V • {(peSg^e^g^t) _ ^ ^ . 

Pw Pg 


(A5) 


5, 


Bu 


Bsh 


+ 


Sg^fi Sg^fi \ 


Ba 


— (1 - </>e) Ps/i + V ■ [(1 - </>e) Vs] = Qh 


(A6) 


Eqn. (A6) is the mass balance relationship for the 


hydrate-soil-composite matrix. Expanding the deriva¬ 


tives in Eqn. (|A6|) and rearranging gives 

(1 - (t>e) 


d 


Psh "t“ Vg * V psh 


Pshgt<Pe 


d 

dt 


Psh 


T Psh^ ’ [(1 Pe) Vg] — Qh 
d 0 

-^ (1 po) ~^Psh Psh~^Pe T Psh^ ’ [(1 Pe) ^s] 

= Qh ■ (A7) 

Using the expression for rate of change of density of 


the composite matrix from Eqn. (27) in Eqn.(A7), we 
get 

■^ 4 >e = - 5 —TTO'- - 5 —T:-P +V- [(l-</>e)Vs]- 

ut ^sh dt Bgji dt Psh 

(A8) 

where, Bgh is bulk modulus of the composite solid. 


Eqn. (A8) describes the rate of change of the effec¬ 


tive porosity due to external stress a and internal fluid 


pore-pressure P. Einally, substituting Eqn. (A8) in Eqn. 


(A5), we obtain 

Sqn . p. Son , 


Bo, 


Bg 


Bn 


Bsh 


dp Id 

dt Bgh dt 


I f d , d 
+ —cr+a — P 

Bgh \dt dt 


Bsh J 

K-k. 


d d 


+ UlIA ) VP 


pw 


Pg 


Eor the hydrate and the soil phases, the mass conser¬ 
vation equations (Eqn. and Eqn. @) described in 
Section 2.1 are used. Adding Eqn. (§ and Eqn. we 
get 

d 

^ [(t>ShPh + (1 - 0) Ps] + V ■ [<i>ShPh + (!-</>) Ps] Vs = qu 


= ( Nh— + ^ - — ) fco As,o Sh (Pe - P) 

Pw Pg Psh , 


Sw,e . Sg^e 
Bw Bg 




a - Pe 


B 


sh 


dt 


Storativity S 


^sh / \ J->g 

T V • \peSw,e P^w,t T V ' [PeSg^e P^g,t 

+ V-Vs = ^ + ^ + ^ 

Pw Pg Psh 


(A9) 




, d „ K ( I 1 , 

1-^ —e-V-— — + —)VP 
Bsh J dt 2 \Pw Pg / 


a X 

h'f 

= ( + MjL-Ph.\ko As,0 Sh (Pe - P). 

Pw Pg Psh J 

c 


(AlO) 


Table 13: Simplified constitutive relationships 


Hydrate reaction 
kinetics 

qg = ko As Mg {Pe - P) 


As — Ag^o X Sh 


J.J Myj . . Mh 

Darcy velocity 

J. kr,cx „ 

^cx,r = —P -VP 

Poe 

k'p^g - k'p^'iD — 0.5 

Effective stress 

a — a' P- aPI 

principle 

isotropic stress a = a' aP 

Linear elastic 

a' = —2Gah e — As/i {tr e) I 

stress-strain law 

€ = 1 (Vu-h V^u) 


isotropic strain —)■ 

^ — V • U = ^Bsh ~Pxsh^ 


a' = -{1/Bm) cr' 


Eqn. (AlO) can be rewritten in a condensed form as 


^ (Pe - P) • (All) 

at at iJjf 

This is the storage equation describing the pressure re¬ 
sponse in a poroelastic hydrate soil. 













































